Regional heatmaps of model performance under extremes
Overview
This notebook will guide you through the code to plot regional heatmaps of model performance (KGE and hit rate) under climatic extremes in R. The following figures will be computed:
Figure 4: KGE regional heatmap for aggregated crops
Supplementary figures 34-36: regional heatmaps of all KGE components for aggregated crops
Supplementary figure 37: regional heatmap of the model hit rate for aggregated crops
Supplementary figure 38: regions used in the heatmap indicated on a world map
The fully processed and filtered data files, available under GGCMI-validation/data/processed/figure_ready_data are required to produce the figures.
We load the necessary libraries
Code
library(tidyverse) # Includes dplyr, ggplot2, readr, etc.library(sjPlot) # to set and configure plotting themelibrary(rnaturalearth) # to create base world map with countrieslibrary(spdep) # spatial functionslibrary(viridis) # color-blind proof color scheme
We configure the plotting theme based on the data visualization blog from Cédric Scherer: https://www.cedricscherer.com/
base_path <-"C:/Users/kobed/Documents/TOAD-PIKlegacy"# Where is the code repo stored?crop_data <-readRDS(file.path(base_path, "extremes_aggr.rds"))
We have to join rnaturalearth regions to create these regional heatmaps
Code
# Load Natural Earth countries datasetworld <-ne_countries(scale ="medium", returnclass ="sf") %>%mutate(ctr = iso_a3_eh) %>%# eh needed to get France and Norway as well filter(name_long !=" Ashmore and Cartier Islands", name_long !="Indian Ocean Territories") %>%# Get rid of duplicate AUS entries dplyr::select(ctr,subregion) %>%st_drop_geometry() %>%# Drop the geometry column because for AUS two different distinct() crop_data <- crop_data %>%mutate(ctr =ifelse(ctr =="JKX", "PAK", ctr)) %>%# Replace "JKX" with "PAK" for Pakistanleft_join(world, by ="ctr") %>%distinct()
Figure 4 - KGE heatmap
First we write a function to calculate the KGE. Note that we work with a customized version (KGE’) such that KGE = 1 - KGE’.
Then we calculate the KGE’ and its components for each climate extreme and model
Code
heatmap_data <- crop_data %>%group_by(subregion, climate_extreme, model) %>%summarise(KGE_components =calculate_kge(difftrend_obs, difftrend_sim) ) %>%mutate(KGE = KGE_components$KGE,r = KGE_components$r,Beta = KGE_components$Beta,Alpha = KGE_components$Alpha ) %>% dplyr::select(-KGE_components) %>%# Remove the list column, no longer neededdistinct()
We define a model order to be visualised
Code
# Model orderorder <-rev(c("benchmark","ensemble","acea", "simplace-lintul5","dssat-pythia","pepic","cygma1p74","ldndc","pdssat","epic-iiasa", "lpj-guess", "promet","lpjml","isam", "crover"))
Finally, we can visualise the heatmap. We use a logarithmic scale to better differentiate low and high values.
Code
heatmap_data$subregion <-factor(heatmap_data$subregion, levels =rev(c("Western Europe", "Eastern Europe", "Southern Europe", "Northern Europe", "Australia and New Zealand", "Seven seas (open ocean)", "Northern America", "Central America", "Carribean", "South America", "Western Asia", "Central Asia", "Southern Asia", "Eastern Asia", "South-Eastern Asia", "Northern Africa", "Western Africa", "Middle Africa", "Eastern Africa", "Southern Africa")) )heatmap_data <- heatmap_data %>%# necessary because we call heatmap data laterna.omit()heatmap_plot <- heatmap_data %>%filter(subregion !="Seven seas (open ocean)", subregion !="Carribean") %>%na.omit() %>%mutate(model =factor(model, levels = order)) %>%# Order the models using the defined orderggplot(aes(x = model, y = subregion, fill =-log(KGE))) +# We have to do - log(KGE) to get colorbar in right ordergeom_tile(color ="white", size =1) +# Add text annotations with the valuesgeom_text(aes(label =round(KGE, 2)), color ="white", size =3) +# Round to 2 decimals# Set labelslabs(x ="", y ="") +scale_fill_viridis(discrete =FALSE, option ="viridis", direction =1, name ="", breaks =c(-max(log(heatmap_data$KGE)), -min(log(heatmap_data$KGE))), # - to get order from low to high performance rightlabels =c("Low Performance", "High Performance") ) +theme(axis.text.x =element_text(angle =45, hjust =1), # Rotate x-axis labels for readabilityaxis.title.x =element_text(size =12),axis.title.y =element_text(size =12),legend.position ="bottom", # Place the legend on the right for better presentationlegend.key.width =unit(5, "cm"), # Adjust the width of the legend keylegend.key.height =unit(1, "cm"), # Adjust the width of the legend keylegend.text =element_text(size =18), # Set the size of the legend textlegend.title =element_text(size =10), # Set the size of the legend title legend.margin =margin(t =10), # Add some margin above the legendlegend.box ="horizontal", # Ensure the legend spreads horizontally,strip.text =element_text(size =18)) +# Increase facet label sizefacet_wrap(~climate_extreme) heatmap_plot
Code
# Save the heatmapggsave(file.path(base_path, "fig5_HeatmapExtremes.png"), heatmap_plot, width =15, height =20, dpi =300)
Supplementary figures 34-36: decomposition of the KGE heatmap
Alpha
On a log-scale
Code
heatmap_plot <- heatmap_data %>%filter(subregion !="Seven seas (open ocean)", subregion !="Carribean") %>%na.omit() %>%mutate(model =factor(model, levels = order)) %>%# Order the models using the defined orderggplot(aes(x = model, y = subregion, fill =log(Alpha))) +geom_tile(color ="white", size =1) +# Add text annotations with the valuesgeom_text(aes(label =round(KGE, 2)), color ="black", size =3) +# Round to 2 decimals# Set labelslabs(x ="", y ="") +scale_fill_gradient2(high ="red",mid ="white",low ="blue",midpoint =0,name ="",limits =c(min(log(heatmap_data$Alpha)), max(log(heatmap_data$Alpha))),breaks =c(min(log(heatmap_data$Alpha)), log(1), max(log(heatmap_data$Alpha))),labels =c("Lower variability", "Similar variability", "Higher variability")) +theme(axis.text.x =element_text(angle =45, hjust =1), # Rotate x-axis labels for readabilityaxis.title.x =element_text(size =12),axis.title.y =element_text(size =12),legend.position ="bottom", # Place the legend on the right for better presentationlegend.key.width =unit(5, "cm"), # Adjust the width of the legend keylegend.key.height =unit(1, "cm"), # Adjust the width of the legend keylegend.text =element_text(size =18), # Set the size of the legend textlegend.title =element_text(size =10), # Set the size of the legend title legend.margin =margin(t =10), # Add some margin above the legendlegend.box ="horizontal", # Ensure the legend spreads horizontally,strip.text =element_text(size =18)) +# Increase facet label sizefacet_wrap(~climate_extreme) heatmap_plot
Code
# Save the heatmapggsave(file.path(base_path, "suppfig34_HeatmapExtremes.png"), heatmap_plot, width =15, height =20, dpi =300)
Beta
Code
heatmap_plot <- heatmap_data %>%filter(subregion !="Seven seas (open ocean)", subregion !="Carribean") %>%na.omit() %>%mutate(model =factor(model, levels = order)) %>%# Order the models using the defined orderggplot(aes(x = model, y = subregion, fill = Beta)) +geom_tile(color ="white", size =1) +# Add text annotations with the valuesgeom_text(aes(label =round(Beta, 2)), color ="black", size =3) +# Round to 2 decimals# Set labelslabs(x ="", y ="") +scale_fill_gradient2(high ="red",mid ="white",low ="blue",midpoint =0,name ="",limits =c(-5 , 5),oob = scales::squish, # squishes out-of-bounds values to limitsbreaks =c(-5, 0, 5), labels =c("Overestimation bias", "No bias", "Underestimation bias")) +theme(axis.text.x =element_text(angle =45, hjust =1), # Rotate x-axis labels for readabilityaxis.title.x =element_text(size =12),axis.title.y =element_text(size =12),legend.position ="bottom", # Place the legend on the right for better presentationlegend.key.width =unit(5, "cm"), # Adjust the width of the legend keylegend.key.height =unit(1, "cm"), # Adjust the width of the legend keylegend.text =element_text(size =18), # Set the size of the legend textlegend.title =element_text(size =10), # Set the size of the legend title legend.margin =margin(t =10), # Add some margin above the legendlegend.box ="horizontal", # Ensure the legend spreads horizontally,strip.text =element_text(size =18)) +# Increase facet label sizefacet_wrap(~climate_extreme) heatmap_plot
Code
# Save the heatmapggsave(file.path(base_path, "suppfig35_HeatmapExtremes.png"), heatmap_plot, width =15, height =20, dpi =300)
r
Code
heatmap_plot <- heatmap_data %>%filter(subregion !="Seven seas (open ocean)", subregion !="Carribean") %>%na.omit() %>%mutate(model =factor(model, levels = order)) %>%# Order the models using the defined orderggplot(aes(x = model, y = subregion, fill = r)) +# We have to do - log(KGE) to get colorbar in right ordergeom_tile(color ="white", size =1) +# Add text annotations with the valuesgeom_text(aes(label =round(r, 2)), color ="black", size =3) +# Round to 2 decimals# Set labelslabs(x ="", y ="") +scale_fill_gradient2(high ="blue",mid ="white",low ="red",midpoint =0,name ="",limits =c(-1 , 1),breaks =c(-1, 0, 1), labels =c("Negatively correlated", "No correlation", "Positively correlated")) +theme(axis.text.x =element_text(angle =45, hjust =1), # Rotate x-axis labels for readabilityaxis.title.x =element_text(size =12),axis.title.y =element_text(size =12),legend.position ="bottom", # Place the legend on the right for better presentationlegend.key.width =unit(5, "cm"), # Adjust the width of the legend keylegend.key.height =unit(1, "cm"), # Adjust the width of the legend keylegend.text =element_text(size =18), # Set the size of the legend textlegend.title =element_text(size =10), # Set the size of the legend title legend.margin =margin(t =10), # Add some margin above the legendlegend.box ="horizontal", # Ensure the legend spreads horizontally,strip.text =element_text(size =18)) +# Increase facet label sizefacet_wrap(~climate_extreme) heatmap_plot
Code
# Save the heatmapggsave(file.path(base_path, "suppfig36_HeatmapExtremes.png"), heatmap_plot, width =15, height =20, dpi =300)
Supplementary figure 37: hit rate heatmap
First we define a function to calculate the hit rate
Code
calculate_hitrate<-function(benchmark, simulated) {# Check where the simulated values are less than 0 hits <- simulated <0# Calculate hit rate as the mean of hits hitrate <-mean(hits) # This will give the proportion of TRUE values (simulated < 0)return(hitrate) # Return all components}
Then we can calculate per region and climate extreme
`summarise()` has grouped output by 'subregion', 'climate_extreme'. You can
override using the `.groups` argument.
Now we can visualise the heatmap
Code
heatmap_data$subregion <-factor(heatmap_data$subregion, levels =rev(c("Western Europe", "Eastern Europe", "Southern Europe", "Northern Europe", "Australia and New Zealand", "Seven seas (open ocean)", "Northern America", "Central America", "Carribean", "South America", "Western Asia", "Central Asia", "Southern Asia", "Eastern Asia", "South-Eastern Asia", "Northern Africa", "Western Africa", "Middle Africa", "Eastern Africa", "Southern Africa")) )heatmap_data <- heatmap_data %>%# necessary because we call heatmap data laterna.omit()heatmap_plot <- heatmap_data %>%filter(subregion !="Seven seas (open ocean)", subregion !="Carribean") %>%na.omit() %>%mutate(model =factor(model, levels = order)) %>%# Order the models using the defined orderggplot(aes(x = model, y = subregion, fill = hitrate)) +# We have to do - log(KGE) to get colorbar in right ordergeom_tile(color ="white", size =1) +# Add text annotations with the valuesgeom_text(aes(label =round(hitrate, 2)), color ="black", size =3) +# Round to 2 decimals# Set labelslabs(x ="", y ="") +scale_fill_viridis(discrete =FALSE, option ="viridis", direction =1, name ="", limits =c(0,1),breaks =c(0, 1),labels =c("None captured", "All captured") ) +theme(axis.text.x =element_text(angle =45, hjust =1), # Rotate x-axis labels for readabilityaxis.title.x =element_text(size =12),axis.title.y =element_text(size =12),legend.position ="bottom", # Place the legend on the right for better presentationlegend.key.width =unit(5, "cm"), # Adjust the width of the legend keylegend.key.height =unit(1, "cm"), # Adjust the width of the legend keylegend.text =element_text(size =18), # Set the size of the legend textlegend.title =element_text(size =10), # Set the size of the legend title legend.margin =margin(t =10), # Add some margin above the legendlegend.box ="horizontal", # Ensure the legend spreads horizontally,strip.text =element_text(size =18)) +# Increase facet label sizefacet_wrap(~climate_extreme) heatmap_plot
Code
# Save the heatmapggsave(file.path(base_path, "suppfig37_HeatmapExtremes.png"), heatmap_plot, width =15, height =20, dpi =300)
Supplementary figure 38: regions
Finally let’s visualise which grid cells belong to which regions from rnaturalearth
Code
world_geometries <-ne_countries(scale ="medium", returnclass ="sf") %>%mutate(ctr = iso_a3_eh) %>%# eh needed to get France and Norway as well filter(name_long !=" Ashmore and Cartier Islands", name_long !="Indian Ocean Territories") %>%# Get rid of duplicate AUS entries dplyr::select(ctr,subregion) %>%distinct() regions_plot <-ggplot(world_geometries) +geom_sf(fill ="gray90", color ="black", size =0.2) +# Base world mapgeom_tile(data = crop_data, aes(x = lon, y = lat, fill = subregion), size =1, alpha =0.6) +scale_fill_viridis_d(option ="C", name ="Subregion") +# Improved color scalelabs(x ="", y ="") +theme(legend.position ="bottom",legend.title =element_text(size =12),legend.text =element_text(size =10),plot.title =element_text(hjust =0.5, size =16) ) regions_plot
Code
# Save the plot ggsave(file.path(base_path, "suppfig38_Regions.png"), regions_plot, width =28, height =15, dpi =300)